Claims 



What is claimed is : 

1. A method of simulating behavior of a flow interacting with an object, the method 
providing a simulated numerical representation in TV dimensions, N > 3, that is 
composed of a plurality of approximated values at a multitude of points in at least 
a part of space where the flow interacts with the object, the approximated values 
being of a physical parameter u of the flow to which is associated a velocity field a 
which determines a preferential direction, by means of a numerical scheme wherein 
at least one spatial p th derivative £> p , p > 1, of the parameter u is approximated 
at the points of the part of space, the method comprising the steps of : 

using for the part of space a discrete N-dimensional grid constructed by N 
families of coordinate lines ; 

computing, in at least one point P of the grid, called the point of computation, 
an approximated value D£ of D p with an error e n , by using values u s of the 
parameter in a collection of grid points, called the stencil 5, and computational 
functions, evaluated with the values u 3 , the computational functions depending 
on the numerical framework in which D p is expressed, 

choosing the computational functions for the approximated value D£ in such 
a way that the approximated value D£ is optimized for the preferential direction, 
and 

wherein the stencil S contains at least one point situated outside all the 
coordinate lines passing through the point of computation P, and the stencil S 
contains at least a first point and a second point, the first point being defined by 
N first coordinate lines of the N families of lines, the second point being defined 
by N second coordinate lines of the TV families of lines, and for at least one family 
Nf of the coordinate lines, the first coordinate line belonging to the family Nf 
is different from and not adjacent to the second coordinate line belonging to the 
same family N/ ; and 

outputting the numerical representation that simulates, for the part of space, 
behavior of the flow interacting with the object. 

2. The method of claim 1, wherein for varying preferential directions, the approxi- 
mation depends continuously on the values u s . 
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3. The method of claim 1, wherein the computing step of D p comprises the steps of : 

providing a local basis B(e u e 2 ,es, . . .) of curvilinear coordinates which has 
the unit vector e\ along the preferential direction, and 

choosing the computational functions so that a contribution to the error e„ 
of at least one pure or one mixed derivative as expressed in the local basis B is 
minimized, while using as a formulation of the values u s of the parameter at each 
of the points of the stencil 5, a truncated Taylor series expansion with respect to 
the point of computation P with an error, called the truncation error e s . 

4. The method of claim 3, wherein the computational functions are individual coef- 
ficients C s and the approximated value D£ is a linear combination of values u s of 
the parameter. 

5. The method of claim 3, wherein the computing step of D£ comprises the steps of : 

using an integral formulation for the computation of the derivative D p with 
computational functions which are fluxes through a control volume, 

computing in at least one volume an approximated value D£ of D p with 
an error e n , where the approximated value is a function of the flux formulation 
employed. 

6. The method of claim 3, wherein the computing step of D£ comprises the steps of : 

subdividing the part of space into elements, containing nodes at which the 
approximation of the physical parameter u is stored, 

defining basis computational functions <£, called interpolation functions, on 
the elements, where the interpolation functions are used to approximate the phys- 
ical value u on the element, 

computing the integral of the derivative D p on the element, using a test 
computational function ip, called weighting function, 

computing in at least one element an approximated value Dp^ of D p with an 
error e n , where the approximated value is a function of the interpolation functions 
<j> and the weighting function ip. 

7. The method of claim 3, wherein the computing step of D£ comprises the steps of : 

subdividing the part of space into elements, containing nodes at which the 
approximation of the physical parameter u is stored, 
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defining basis computational functions 0, called interpolation functions, on 
the elements, where the interpolation functions axe used to approximate the phys- 
ical value u on the element, 

computing the integral I e i of the derivative D p on volumes, and distributing 
parts ctilei to nodes z, where a, represent computational functions called distribu- 
tion coefficients. 

8. The method of claim 3, wherein the computing step of D£ comprises the steps of : 

subdividing the part of space into elements, containing nodes at which the 
approximation of the physical parameter u is stored, 

defining basis computational functions called interpolation functions, on 
the elements, where the interpolation functions are used to approximate the phys- 
ical value u on the element, 

computing fluxes / at the surfaces of volumes, and distributing parts aif to 
nodes i, where represent computational functions called distribution coefficients, 

9. The method of claim 1, wherein the computing step of D£ comprises the steps of : 

using a representation of the numerical solutions in Fourier components, and 

choosing in the approximated value D£ the computational functions in such 
a way that the Fourier components are optimized for certain directions which are 
related to the velocity a, while using the values u s of the parameter at each of the 
points of the stencil S in the Fourier representation. 

10. The method of claim 9, wherein the computational functions are individual coef- 
ficients C 3 and the approximated value D£ is a linear combination of values u s of 
the parameter. 

11. The method of claim 9, wherein the computing step of comprises the steps of : 

using an integral formulation for the computation of the derivative D p with 
computational functions which axe fluxes through a control volume, 

computing in at least one volume an approximated value of D p with 
an error e n , where the approximated value is a function of the flux formulation 
employed. 

12. The method of claim 9, wherein the computing step of D£ comprises the steps of : 

subdividing the part of space into elements, containing nodes at which the 
approximation of the physical parameter u is stored, 
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defining basis computational functions <£, called interpolation functions, on 
the elements, where the interpolation functions are used to approximate the phys- 
ical value u on the element, 

computing the integral of the derivative D p on the element, using a test 
computational function -0, called weighting function, 

computing in at least one element an approximated value Dp* of D p with an 
error c n , where the approximated value is a function of the interpolation functions 
<j> and the weighting function tp. 

13. The method of claim 9, wherein the computing step of D£ comprises the steps of : 

subdividing the part of space into elements, containing nodes at which the 
approximation of the physical parameter u is stored, 

defining basis computational functions 0, called interpolation functions, on 
the elements, where the interpolation functions are used to approximate the phys- 
ical value u on the element, 

computing the integral I et of the derivative D p on volumes, and distributing 
parts OiiI e i to nodes i, where a* represent computational functions called distribu- 
tion coefficients. 

14. The method of claim 9, wherein the computing step of comprises the steps of : 

subdividing the part of space into elements, containing nodes at which the 
approximation of the physical parameter u is stored, 

defining basis computational functions 0, called interpolation functions, on 
the elements, where the interpolation functions are used to approximate the phys- 
ical value u on the element, 

computing fluxes / at the surfaces of volumes, and distributing parts a»/ to 
nodes z, where a* represent computational functions called distribution coefficients, 

15. The method of claim 1, wherein the TV-dimensional grid is expressed in a coordi- 
nate system which is chosen from the group consisting of : rectangular coordinates, 
spherical coordinates, cylindrical coordinates, parabolic cylindrical coordinates, 
paraboloidal coordinates, elliptic cylindrical coordinates, prolate spheroidal coor- 
dinates, oblate spheroidal coordinates, bipolar coordinates, toroidal coordinates, 
conical coordinates, confocal ellipsoidal coordinates and confocal paraboloidal co- 
ordinates. 

16. The method of claim 1, wherein the TV-dimensional grid is chosen from the group 
consisting of : a grid with non-uniform mesh spacing, a grid which is moving, a 
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grid which is deforming, a grid which is rotating, and a grid which is staggered, 
and any combination thereof. 

17. The method of claim 1, wherein the spatial p th derivative D p is a pure derivative 

18. The method of claim 1, wherein p = 1. 

19. The method of claim 1, wherein p = 2. 

20. The method of claim 1, wherein the spatial p th derivative D p is a mixed derivative 
^Idf^dg? . . . with Pl +p 2 + ...=p. 

21. The method of claim 3, computing an approximated value of D x , where the 
approximated value is denoted by 

{ (Ux)i,j,k,...\ 

( u v)i,j,k,... 



V 



(27) 



wherein 



+ a n _iU i+n _i j:?)A:v .. 4- a n n i+n)J)iki ... + Tx j , 

(^y)*,^,... = | a -m ( u i~mj~m,k,... ~ Ui-mJ.k,...) 

H- a_ m+ i(^i_ m+ i }:? _ m _j-i jA ; v .. — l/i-m+i^A:,...) + • • - 
H-an-l(Wt+n-l 5 j+n-l,A:,... ~ ,j,A;,...) 
+ a n (u^ fn ,j+n,A 1 ... - Wi+nj,*,...) + | , 

C^z)*,.?,/:,... = I ^-m^j-mj-mjfc-m,,.. ~ ^t-mj-m,*,...) 

+ an-l^t+n-lj-fn-l^+n-l,... ~ ^t+n-lj+n-l^,...) 

4- a„(^ +nj j +n ^ +nv .. - w i+n)j+n>fcj ...) +T*}, 



(28) 



wherein a_ m 7^ 0, a n ^ 0, m and n are given integers, m 4- n > 0, and m 4- n > 
1 if m * n = 0, the terms Trc, Ty, Tz, . . . represent the degrees of freedom 
which are used in the optimization of the approximated value JDj 4 , and where 
indices (i, j, /c, . . .) define the point of computation P on the iV-dimensional grid, 
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and Ax, Ay, Az, . . . denote the mesh spacings of the JV-dimensional grid in each 
coordinate direction. 

22. The method of claim 21, computing an approximated value of D\ = ^ u /de\ 
by a discretization of order M, in which the terms d M ^ l u/de^ 2 de^ . . . with 
M 2 + M 3 4- . . . = M -h 1 are eliminated in the approximated value D*. 

23. The method of claim 21, computing an approximated value of Di = ® u /de\ 
by a discretization of order M, in which the terms d M ~ { ' 1 u/de^ fl de^ r2 de^ . . . with 
Mi + Mi + M 3 + . . . = M + 1 and Mi < M + 1 are eliminated in the approximated 
value in the case that e\ is along the x-axis or along diagonals. 

24. The method of claim 1, computing an approximated value of Z>i, where the 
approximated value is denoted by 

V : 



(29) 



wherein 

-fa n _iUi +n _i j:7j A;,... + a n u i+nijik ^ I , 

( U y)id,k,... = | fl -m(^t-m,j-m ) fc ) ... ~~ ^t-mj,/:,...) 

+ fl-m+l(^i-m+lj-m+l,A:,... ~~ ^i-m+lj ,*,...) + • • • 

-j- a n _i(-Ui +n _i }J+n _i^,... — Wt+n-ij,*,...) 

(^2 )*,.?,/:,... = "^y^ ^ a — m(^»— m,j— m,fc— m,... ^i—mj— m,A:,...) 

+ a_ m+ i('Ui_ m+ i ) j_ m+ i j A : - m+ i,... — Ui-m+ij-m+i^,...) + . . . 
+ On-l^i+n-lj+n-l^+n-l,... Ui+n-lJ+n-l,*,...) 

: (so) 

wherein a_ m ^ 0, a n ^ 0, m and n are given integers, ra + n > 0, and ra + n > 1 
if m * n = 0, and where indices (i, j, A:, . . .) define the point of computation P 
on the TV-dimensional grid, and Ax, Ay, Az, . . . denote the mesh spacings of the 
iV-dimensional grid in each coordinate direction. 
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25. The method of claim 21, wherein 

(«x)ij,fc,... = 2^ ~ U i-U,k,..) + Tx} , 

( u v)ij,k, -- = { \( u i+l J+l.*,... ~ Wi+l,i,A.... + «t-l j,*,... - «t-l j-l,fc,...) + ^} > 

(«r)tj,fc,... = ^ jg^+lj+Lfc+l, .. - "i+lj+l,*,... 
5 + Ui-ij-!,*,... - tii-ij-Lfc-i,...) + Tz} , 

5 (31) 

26. The method of claim 21, wherein 

- 8^-1^-1,*,... + 8^-1^,... 
+ Siij+ij+i,*,... - Sut+ij,*,... 

- 1^+2^+2,*,... + ^+2^,*,. J + Tj/j , 
{u z )i^ K „ = — | — (U i _ 2j j-2,fc-2 s ... - Wt-2J-2,k,... 

+ 8^1+1^+1^+1,... — 8^+1^+1^,., 



— Ui+2j+2,k+2 t ... + ^+2,^+2,*:,...) + Tz j , 



20 27. The method of claim 21, wherein 

+ |(tt<-2J-2 I Jb,...-tti-2J 1 * l ...)+Ty^ , 
25 + ^K-2,j-2,fc-2,... - ^-2^-2,*,..) + TZ^ , 



(32) 



(33) 
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28. A simulation method according to claim 1, in three dimensions, obtaining the 
approximation with order M of the derivative D v = ^^/de^de^de? 3 wit ^ 
Pi +P2 +P3 = P on a S r ^d °f given extent from the output of the program generate- 
discretizations which is given in appendices 1-5. 

5 29. A simulation method according to claim 28, wherein the input parameter optimize 
has the value 1 or 2. 

30. A simulation method according to claim 28, wherein the input parameter order 
has the value 1. 

31. The method of claim 3, computing in three dimensions an approximation of the 
10 second derivative D 2 = ^ u /de\y where the approximation D$ is expressed in the 

terms u xx , Uyy, u zz , u xy , u yz and u zx which are given by 

(V>xx)i,j,k = J.* 2u ij>* + U i~hj,k + TxX ) > 

(«w)w.A = Kj+i,* - 2ui dik + u id - hk + Tyy) , 
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-Ui- lij)k - n^-i^ -h Ui- ld - 1)k ) + Txy} , 

( u yz)iJ,k = AyA* - Wi-fi^+i - + U i+ i J)A: + 

~ u ij,k-l + UtJ-l,*-l + 

20 - Wi-ij-i,* - t/i-i^-i 4- Ui- Ujtk ) + T$/2 j , 

— Ui-^/c - U iijik -l + ^-1^-1 + 

Ui-ij- lik -i - Ui- hj - hk - iiij-i.fc-1 + Wtj-i,*) + T2X j (34) 

25 wherein the terms Txx, Txy, Txz, Tyy, Tyz and Tzz represent the degrees of 

freedom which are used in the optimization of the approximated value D£> and 
where indices (2, j, k) define the point of computation P on the three-dimensional 
grid, and Ax, Ay, Az denote the mesh spacings of the three-dimensional grid in 
each coordinate direction. 

30 32. The method of claim 1, wherein D£ = £ n L n D^ n where L n are constants and each 
£>p n is a function of values u 3 of the parameter in a collection of grid points, called 
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stencil S n , with individual computation functions, which depend on the numerical 
framework in which D£ n is expressed, and wherein in the approximation D£ n i tne 
computation functions are chosen in such a way that the approximated value D£ l n 
is optimized for the preferential direction. 

33. The method of claim 1, wherein D£ = J2 n L nD^ n where L n are limiting functions 
of the values u 3 of the stencil S, and at least one D£ l n is a function of values u 3 
of the parameter in a collection of grid points, called stencil S ni with individual 
computation functions, which depend on the numerical framework in which D£ n 
is expressed, and wherein in the approximation D£ n , the computation functions 
are chosen in such a way that said approximated value D£ n is optimized for the 
preferential direction. 

34. The method of claim 1, wherein the stencil S is chosen from the group consisting 
of : upwind discretization stencils, centered discretization stencils, and discretiza- 
tion stencils which are a combination of at least one upwind discretization stencil 
and at least one centered discretization stencil. 

35. The method of claim 1, wherein the numerical discretization is a non-linear dis- 
cretization. 

36. The method of claim 1, wherein the numerical scheme is chosen from the group 
consisting of : the Lax-Wendroff scheme, the Lax-Priedrich scheme, the MacCor- 
mack scheme, the leap-frog scheme, the Crank-Nicholson scheme, the Stone-Brian 
scheme, the box scheme, Henn's scheme, the QUICK scheme, the k, scheme, the 
Flux Corrected Transport (FTC) scheme, the family of ENO schemes, schemes 
in the class of the Piecewise Parabolic Method (PPM), multi-level schemes, and 
schemes obtained with the fractional step method and variants thereof. 

37. The method of claim 1, wherein the numerical scheme includes the discretization 
of a plurality of equations. 

38. The method of claim 1 for the numerical simulation of physical phenomena which 
are modeled by the Navier-Stokes equations with equation (s) of state. 

39. The method of claim 1 for the numerical simulation of physical phenomena which 
are modeled by the Euler equations with equation(s) of state. 

40.. The method of claim 1 for the numerical simulation of physical phenomena which 
are modeled by the magneto-hydrodynamic equations with equation(s) of state. 
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41. The method of claim 1 in combination with at least one model for the physical 
phenomena chosen from the group consisting of : equations to model turbulence, 
equations to model chemical reactions, equations to model electromagnetic phe- 
nomena, equations to model multiphase flow and equations to model multi-physics 

5 phenomena. 

42. The method of claim 1 for the numerical simulation of physical phenomena in com- 
bination with at least one acceleration technique chosen from the group consisting 
of : local time-stepping, multi-grid, GMRES and preconditioning. 

43. The method of claim 1 for the simulation of a physical phenomenon which includes 
10 a material flow. 

44. The method of claim 1 for the simulation of a physical phenomenon which includes 
a material object interacting with material flow. 

45. The method of claim 1 for the simulation of a physical phenomenon which includes 
a vehicle. 

15 46. The method of claim 1 for the simulation of a physical phenomenon which includes 
a rotating blade. 

47. The method of claim 1 for the simulation of at least a part of the atmosphere of 
the earth. 

48. The method of claim 1 for the numerical simulation of oil recovery. 

. 20 49. The method of claim 1 for the numerical simulation of physical phenomena which 
include combustion. 

50. A data processing system programmed to implement a simulation method accord- 
ing to claim 1. 

51. A computer program that can be loaded in a data processing system so as to 
25 implement a method according to claim 1. 

52. A digital storage computer-readable medium containing a stored computer pro- 
gram that is configured to be loadable in a data processing system so as to imple- 
ment a method according to claim 1. 
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